Frog Phenology Statisics V3 Mar24
Xav Harrison
2024-02-20
#Libraries and Data
##Libraries
library(brms)## Loading required package: Rcpp## Loading 'brms' package (version 2.21.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:stats':
##
## ar library(factoextra)## Loading required package: ggplot2## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa library(FactoMineR)
library(corrplot)## corrplot 0.92 loaded library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors library(cowplot)##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stamp library(naniar)
library(rstatix)##
## Attaching package: 'rstatix'
##
## The following object is masked from 'package:stats':
##
## filter library(RColorBrewer)
library(ggeffects)##
## Attaching package: 'ggeffects'
##
## The following object is masked from 'package:cowplot':
##
## get_title plotopts<-theme_bw() +theme(axis.text=element_text(size=12),axis.title = element_text(size=15),legend.text = element_text(size=15),legend.title = element_text(size=15),strip.text.x = element_text(size=16))Data
#Pond Data
pond_wide<-read.csv('Pond Wide Assembled Data.csv',header=T)
dim(pond_wide)## [1] 3944 58#Land Cover Data
land_cover_wide<-read.csv('Land Cover Wide Assembled Data.csv',header=T,row.names = 1)
#Spawn Climate
spawn_climate<-read.csv('Spawn Climate Assembled Data.csv',header=T,row.names = 1)
#Winter Climate
winter_climate<-read.csv('Winter Climate Assembled Data.csv',header=T,row.names = 1)
#Tadpole Climate
tadpole_climate<-read.csv('Tadpole Climate Assembled Data.csv',header=T,row.names = 1)
## Year Coverage
pond_wide %>% group_by(SITECODE) %>% reframe(min=min(survey_year),max=max(survey_year))## # A tibble: 10 × 3
## SITECODE min max
## <chr> <int> <int>
## 1 T01 1994 2012
## 2 T02 1994 2014
## 3 T03 1994 2011
## 4 T04 1994 2015
## 5 T05 1994 2015
## 6 T06 1994 2010
## 7 T08 1994 1994
## 8 T09 1995 2001
## 9 T10 2001 2015
## 10 T11 1995 2012Housekeeping
#Variable Labels
label_df<-data.frame(old=c("year","arable_area","grassland_area","tmax","urban_area","spawn_no3n","forest_area","tmin","af","freshwater_area","other_area","spawn_nh4n","rain"),new=c("Year","Arable Area","Grassland Area","Maximum Temperature","Urban Area","Spawn Nitrate","Forest Area","Minimum Temperature","Air Frost","Freshwater Area","Other Area","Spawn Ammonium","Rain"))Dataset Assembly
SPAWNDATE data frame (Winter Climate)
#Creating a data frame containing the earliest spawning date, the highest ammonium concentration breeding adults were exposed to???, and the highest nitrate concentration breeding adults were exposed to??? in each ECN location, pond, and year
pond_wide %>% group_by(SITECODE,LCODE,survey_year) %>% summarize(spawn=mean(spawn_yearday,na.rm=TRUE),
prespawn_nh4n=mean(nh4n_numeric[survey_yearday<min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE),
prespawn_no3n=mean(no3n_numeric[survey_yearday<min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE)
) -> spawndate_df## Warning: There were 78 warnings in `summarize()`.
## The first warning was:
## ℹ In argument: `prespawn_nh4n = mean(...)`.
## ℹ In group 18: `SITECODE = "T01"`, `LCODE = 1`, `survey_year = 2011`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 77 remaining warnings.## `summarise()` has grouped output by 'SITECODE', 'LCODE'. You can override using
## the `.groups` argument. colnames(spawndate_df)[colnames(spawndate_df) == "survey_year"] <- "year"
colnames(spawndate_df)[colnames(spawndate_df) == "SITECODE"] <- "location"
#Adding land cover data
spawndate_df = merge(x=spawndate_df,y=land_cover_wide,by=c("year", "location"),all.x=TRUE)
#Adding winter climate data
spawndate_df = merge(x=spawndate_df,y=winter_climate,by=c("year", "location"),all.x=TRUE)
#Replacing infinity values with NA
spawndate_df %>% replace_with_na_all(condition = ~.x == -Inf)-> spawndate_df
spawndate_df %>% replace_with_na_all(condition = ~.x == Inf)-> spawndate_dfHATCHDATE data frame (Spawn Climate)
#Creating a data frame containing the earliest hatching date, the highest ammonium concentration spawn were exposed to???, and the highest nitrate concentration spawn were exposed to??? in each ECN location, pond, and year
pond_wide %>% group_by(SITECODE,LCODE,survey_year) %>% summarize(hatch=mean(hatch_yearday,na.rm=TRUE),
spawn_nh4n=mean(nh4n_numeric[survey_yearday<min(hatch_yearday,na.rm=TRUE) & survey_yearday>min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE),
spawn_no3n=mean(no3n_numeric[survey_yearday<min(hatch_yearday,na.rm=TRUE) & survey_yearday>min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE)
)-> hatchdate_df## Warning: There were 222 warnings in `summarize()`.
## The first warning was:
## ℹ In argument: `spawn_nh4n = mean(...)`.
## ℹ In group 1: `SITECODE = "T01"`, `LCODE = 1`, `survey_year = 1994`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 221 remaining warnings.## `summarise()` has grouped output by 'SITECODE', 'LCODE'. You can override using
## the `.groups` argument. colnames(hatchdate_df)[colnames(hatchdate_df) == "survey_year"] <- "year"
colnames(hatchdate_df)[colnames(hatchdate_df) == "SITECODE"] <- "location"
#Adding land cover data
hatchdate_df = merge(x=hatchdate_df,y=land_cover_wide,by=c("year", "location"),all.x=TRUE)
#Adding data on the climate spawn were exposed to
hatchdate_df = merge(x=hatchdate_df,y=spawn_climate,by=c("year", "location"),all.x=TRUE)
#Replacing infinity values with NA
hatchdate_df %>% replace_with_na_all(condition = ~.x == -Inf)-> hatchdate_df
hatchdate_df %>% replace_with_na_all(condition = ~.x == Inf)-> hatchdate_dfPERCDEAD data frame
#Creating a data frame containing the highest percentage of spawn found dead, the highest ammonium concentration spawn were exposed to???, and the highest nitrate concentration spawn were exposed to??? in each ECN location, pond, and year
pond_wide %>% group_by(SITECODE,LCODE,survey_year) %>% dplyr::summarize(percdead=mean(percdead_numeric, na.rm=TRUE),
spawn_nh4n=mean(nh4n_numeric[survey_yearday<min(hatch_yearday,na.rm=TRUE) & survey_yearday>min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE),
spawn_no3n=mean(no3n_numeric[survey_yearday<min(hatch_yearday,na.rm=TRUE) & survey_yearday>min(spawn_yearday,na.rm=TRUE)], na.rm=TRUE)
)-> percdead_df## Warning: There were 222 warnings in `dplyr::summarize()`.
## The first warning was:
## ℹ In argument: `spawn_nh4n = mean(...)`.
## ℹ In group 1: `SITECODE = "T01"`, `LCODE = 1`, `survey_year = 1994`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 221 remaining warnings.## `summarise()` has grouped output by 'SITECODE', 'LCODE'. You can override using
## the `.groups` argument. colnames(percdead_df)[colnames(percdead_df) == "survey_year"] <- "year"
colnames(percdead_df)[colnames(percdead_df) == "SITECODE"] <- "location"
#Adding land cover data
percdead_df = merge(x=percdead_df,y=land_cover_wide,by=c("year", "location"),all.x=TRUE)
#Adding data on the climate spawn were exposed to - should I do winter climate instead/as well, a measure of adult fitness?
percdead_df = merge(x=percdead_df,y=spawn_climate,by=c("year", "location"),all.x=TRUE)
#Replacing infinity values with NA
percdead_df %>% replace_with_na_all(condition = ~.x == -Inf)-> percdead_df
percdead_df %>% replace_with_na_all(condition = ~.x == Inf)-> percdead_df
head(percdead_df)## # A tibble: 6 × 17
## year location LCODE percdead spawn_nh4n spawn_no3n forest_area arable_area
## <int> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1994 T01 1 32.9 NaN NaN 608706. 9789443.
## 2 1994 T02 1 15 NaN NaN 1840721. 71881.
## 3 1994 T02 2 NaN NaN NaN 1840721. 71881.
## 4 1994 T03 1 0.167 NaN NaN 1808634. 3341825.
## 5 1994 T04 4 1.43 NaN NaN 22843. 0
## 6 1994 T04 1 0 NaN NaN 22843. 0
## # ℹ 9 more variables: grassland_area <dbl>, freshwater_area <dbl>,
## # urban_area <dbl>, other_area <dbl>, tmax <dbl>, tmin <dbl>, af <dbl>,
## # rain <dbl>, sun <dbl> nrow(percdead_df)## [1] 224 #
ggplot(percdead_df, aes(y=percdead, x=year)) + geom_point(shape=21,size=5,aes(fill=factor(LCODE)))## Warning: Removed 27 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot(percdead_df, aes(y=percdead, x=spawn_no3n)) + geom_point(shape=21,size=5,aes(fill=factor(LCODE)))## Warning: Removed 189 rows containing missing values or values outside the scale range
## (`geom_point()`).
# percdead_yearSample Sizes and Filtering
#HATCH
hatchdate_df %>% group_by(location) %>%reframe(hatch=length(hatch[which(!is.na(hatch))]))## # A tibble: 10 × 2
## location hatch
## <chr> <int>
## 1 T01 7
## 2 T02 23
## 3 T03 13
## 4 T04 45
## 5 T05 21
## 6 T06 10
## 7 T08 1
## 8 T09 4
## 9 T10 2
## 10 T11 26 hatchdate_df<-filter(hatchdate_df,location %in% c("T01","T02","T03","T04","T05","T06","T11") & hatch>0)
nrow(hatchdate_df)## [1] 145 table(hatchdate_df$location)##
## T01 T02 T03 T04 T05 T06 T11
## 7 23 13 45 21 10 26 hatchdate_df$location<-factor(hatchdate_df$location)
#SPAWN
spawndate_df %>% group_by(location) %>%reframe(spawn=length(spawn[which(!is.na(spawn))]))## # A tibble: 10 × 2
## location spawn
## <chr> <int>
## 1 T01 17
## 2 T02 23
## 3 T03 19
## 4 T04 52
## 5 T05 22
## 6 T06 13
## 7 T08 1
## 8 T09 5
## 9 T10 4
## 10 T11 29 spawndate_df<-filter(spawndate_df,location %in% c("T01","T02","T03","T04","T05","T06","T11") & spawn>0)
nrow(spawndate_df)## [1] 175 spawndate_df$location<-factor(spawndate_df$location)
#PERCDEAD
percdead_df %>% group_by(location) %>%reframe(percdead=length(percdead[which(!is.na(percdead))]))## # A tibble: 10 × 2
## location percdead
## <chr> <int>
## 1 T01 16
## 2 T02 22
## 3 T03 14
## 4 T04 67
## 5 T05 18
## 6 T06 13
## 7 T08 0
## 8 T09 5
## 9 T10 14
## 10 T11 28 percdead_df<-filter(percdead_df,location %in% c("T01","T02","T03","T04","T05","T06","T11") & percdead>=0)
nrow(percdead_df)## [1] 178 table(percdead_df$location)##
## T01 T02 T03 T04 T05 T06 T11
## 16 22 14 67 18 13 28 percdead_df$location<-factor(percdead_df$location)PCA and MODELS
SPAWNDATE
Mean Trend Over Time
#Sample Size
spawndate_df %>% filter(!is.na(spawn)) %>% nrow()## [1] 175spawnmeans<- spawndate_df %>% group_by(year,location) %>% reframe(spawn=mean(spawn,na.rm=T))
spawn_siteplot1<-ggplot(spawnmeans,aes(x=year,y=spawn))+ geom_smooth(method="lm") + geom_point(shape=21,size=5,aes(fill=location)) + facet_wrap(.~location,nrow=4) + guides(fill="none") + labs(x="Year",y="Spawn Date (Days 1st Jan)") + plotopts + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2")
spawn_siteplot1## `geom_smooth()` using formula = 'y ~ x'
# Correlations
spawndate_df %>% group_by(location) %>% cor_test(spawn,year,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))## # A tibble: 7 × 10
## location var1 var2 cor statistic p conf.low conf.high method p.adj
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 T01 spawn year -0.7 -3.76 0.0019 -0.882 -0.324 Pearson 0.0133
## 2 T02 spawn year -0.34 -1.63 0.118 -0.657 0.0891 Pearson 0.275
## 3 T03 spawn year -0.068 -0.283 0.781 -0.507 0.398 Pearson 0.781
## 4 T04 spawn year -0.26 -1.92 0.0605 -0.499 0.0116 Pearson 0.212
## 5 T05 spawn year 0.28 1.28 0.215 -0.165 0.625 Pearson 0.301
## 6 T06 spawn year 0.16 0.554 0.59 -0.425 0.656 Pearson 0.688
## 7 T11 spawn year -0.25 -1.32 0.199 -0.562 0.133 Pearson 0.301#LMS
#Models
spawn_mods<-spawndate_df %>% group_by(location) %>% do(tidy(lm(spawn~year,data=.))) %>% filter(term=="year") #%>% mutate(p=p.value) %>% mutate(p.adj=p.adjust(p,method="BH"))
spawn_mods$p.adj<-p.adjust(p=spawn_mods$p.value,"BH")
#Annotations
ann_text <- data.frame(year = 2005,spawn = 100,lab = "cor = -0.7 \np.adj = 0.01",
location = "T01")
spawn_siteplot2<-spawn_siteplot1 + geom_text(data = ann_text,label = ann_text$lab,size=6)
spawn_siteplot2## `geom_smooth()` using formula = 'y ~ x'
#Models
spawn_site_model1<-glm(spawn ~ location*year,data=spawndate_df)
summary(spawn_site_model1)##
## Call:
## glm(formula = spawn ~ location * year, data = spawndate_df)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3029.3824 1053.7266 2.875 0.00459 **
## locationT02 -1656.2168 1263.0377 -1.311 0.19162
## locationT03 -2697.8938 1397.4266 -1.931 0.05529 .
## locationT04 -1966.8514 1148.6924 -1.712 0.08878 .
## locationT05 -3765.8058 1274.0553 -2.956 0.00359 **
## locationT06 -3561.9439 1563.8236 -2.278 0.02406 *
## locationT11 -1909.4359 1355.9569 -1.408 0.16100
## year -1.4853 0.5263 -2.822 0.00538 **
## locationT02:year 0.8383 0.6308 1.329 0.18575
## locationT03:year 1.3530 0.6980 1.938 0.05433 .
## locationT04:year 1.0042 0.5737 1.751 0.08194 .
## locationT05:year 1.8642 0.6361 2.930 0.00388 **
## locationT06:year 1.7865 0.7814 2.286 0.02354 *
## locationT11:year 0.9582 0.6769 1.415 0.15886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 113.0278)
##
## Null deviance: 117574 on 174 degrees of freedom
## Residual deviance: 18197 on 161 degrees of freedom
## AIC: 1339.4
##
## Number of Fisher Scoring iterations: 2 plot(ggeffects::ggeffect(spawn_site_model1,term= ~year*location)) 
confint(spawn_site_model1)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 9.641161e+02 5094.6485864
## locationT02 -4.131725e+03 819.2916000
## locationT03 -5.436800e+03 41.0119295
## locationT04 -4.218247e+03 284.5443839
## locationT05 -6.262908e+03 -1268.7034183
## locationT06 -6.626982e+03 -496.9058782
## locationT11 -4.567063e+03 748.1907367
## year -2.516893e+00 -0.4536957
## locationT02:year -3.980331e-01 2.0745529
## locationT03:year -1.507074e-02 2.7210279
## locationT04:year -1.201660e-01 2.1286619
## locationT05:year 6.173680e-01 3.1109842
## locationT06:year 2.550412e-01 3.3179075
## locationT11:year -3.685586e-01 2.2848687Mean Trend with AF Days
spawn_af_plot1<-ggplot(spawndate_df,aes(x=af,y=spawn))+ geom_smooth(method="lm") + geom_point(shape=21,size=5,aes(fill=location)) + facet_wrap(.~location,nrow=4) + guides(fill="none") + labs(x="Air Frost Days",y="Spawn Date (Days 1st Jan)") + plotopts + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2")
spawn_af_plot1## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Correlations
spawndate_df %>% group_by(location) %>% cor_test(spawn,af,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))## # A tibble: 7 × 10
## location var1 var2 cor statistic p conf.low conf.high method p.adj
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 T01 spawn af -0.0041 -0.0159 0.987 -0.484 0.477 Pearson 0.987
## 2 T02 spawn af 0.27 1.27 0.217 -0.162 0.612 Pearson 0.401
## 3 T03 spawn af 0.23 0.994 0.334 -0.246 0.622 Pearson 0.468
## 4 T04 spawn af 0.25 1.82 0.0754 -0.0260 0.488 Pearson 0.321
## 5 T05 spawn af 0.4 1.79 0.0917 -0.0689 0.722 Pearson 0.321
## 6 T06 spawn af 0.36 1.27 0.229 -0.240 0.759 Pearson 0.401
## 7 T11 spawn af 0.061 0.320 0.751 -0.312 0.419 Pearson 0.876Mean Trend with Temperature
spawn_tmax_plot1<-ggplot(spawndate_df)+ geom_smooth(method="lm",aes(x=tmax,y=spawn,group=location)) + geom_point(shape=21,size=5,aes(fill=location,x=tmax,y=spawn)) + guides(fill="none") + labs(x="Mean Winter Temperature",y="Spawn Date (Days 1st Jan)") + plotopts + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2")
spawn_tmax_plot1## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Correlations
spawndate_df %>% group_by(location) %>% cor_test(spawn,tmax,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))## # A tibble: 7 × 10
## location var1 var2 cor statistic p conf.low conf.high method p.adj
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 T01 spawn tmax -0.31 -1.28 0.219 -0.691 0.196 Pearson 0.307
## 2 T02 spawn tmax -0.45 -2.33 0.0298 -0.729 -0.0508 Pearson 0.0695
## 3 T03 spawn tmax -0.19 -0.804 0.433 -0.594 0.288 Pearson 0.449
## 4 T04 spawn tmax -0.34 -2.56 0.0135 -0.561 -0.0746 Pearson 0.0695
## 5 T05 spawn tmax -0.52 -2.52 0.0218 -0.789 -0.0892 Pearson 0.0695
## 6 T06 spawn tmax -0.37 -1.31 0.217 -0.764 0.230 Pearson 0.307
## 7 T11 spawn tmax -0.15 -0.768 0.449 -0.487 0.233 Pearson 0.449PCA
#Removing sun
spawndate_nosun <- spawndate_df %>% dplyr::select(-sun)
#Removing all rows containing NA values
spawndate_nona <- na.omit(spawndate_nosun)
#Removing response variable
spawndate_noresp <- spawndate_nona %>% dplyr::select(-spawn)
#Removing location variables
spawndate_for_pca <- spawndate_noresp %>% dplyr::select(-location)
spawndate_for_pca <- spawndate_for_pca %>% dplyr::select(-LCODE)
head(spawndate_for_pca) ## # A tibble: 6 × 13
## year prespawn_nh4n prespawn_no3n forest_area arable_area grassland_area
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1994 0 0.0145 1492094. 5194537. 12178451.
## 2 1996 0.25 0.051 21442. 0 19291786.
## 3 1996 0.04 0.044 21442. 0 19291786.
## 4 1996 0.08 0.0625 21442. 0 19291786.
## 5 1997 0.05 0.02 20741. 0 19292486.
## 6 1997 0.25 0.039 20741. 0 19292486.
## # ℹ 7 more variables: freshwater_area <dbl>, urban_area <dbl>,
## # other_area <dbl>, tmax <dbl>, tmin <dbl>, af <dbl>, rain <dbl>#Variable Rename
spawndatepcalabs<-label_df$new[match(colnames(spawndate_for_pca),label_df$old)]
spawndatepcalabs[c(2,3)]<-c("Pre-Spawn Ammonium","Pre-Spawn Nitrate")
colnames(spawndate_for_pca)<-spawndatepcalabs
#Performing Principal Components Analysis
spawndate_res_pca <- PCA(spawndate_for_pca, graph = FALSE)
#Contributions
dimdesc(spawndate_res_pca, axes = 1:2)## $Dim.1
##
## Link between the variable and the continuous variables (R-square)
## =================================================================================
## correlation p.value
## Forest Area 0.9629334 1.997195e-13
## Arable Area 0.9577990 7.623573e-13
## Freshwater Area 0.9219618 4.132568e-10
## Urban Area 0.9191159 5.943078e-10
## Other Area 0.8069561 3.284910e-06
## Minimum Temperature 0.6803874 3.535380e-04
## Maximum Temperature 0.6407077 9.885594e-04
## Air Frost -0.5568433 5.782057e-03
## Grassland Area -0.9579159 7.408580e-13
##
## $Dim.2
##
## Link between the variable and the continuous variables (R-square)
## =================================================================================
## correlation p.value
## Air Frost 0.6152854 0.0017786579
## Minimum Temperature -0.5681398 0.0046807595
## Maximum Temperature -0.6429724 0.0009358000
## Rain -0.6650735 0.0005351047#Producing a rotation plot
spawndate_rotationplot<-fviz_pca_var(spawndate_res_pca, col.var="black",
repel = TRUE) +
labs(x = "PC1", y = "PC2", title = "") + guides(color="none")Spawndate Data
#Creating a data frame containing the co-ordinates of PC1 and PC2
spawndate_res_pca$ind$coord## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 8.26981844 -1.8260942 1.39323957 1.52336216 -1.11175053
## 2 -2.00017040 1.9034884 2.43094661 1.10103625 -0.82950343
## 3 -1.53808478 1.2877144 0.86106097 -0.04125161 -0.18623237
## 4 -1.60423205 1.5265109 1.39354157 0.01797176 0.06829241
## 5 -1.16734789 0.4468061 0.42354312 0.06710644 -0.56053220
## 6 -1.59171770 1.1205533 2.08639707 1.04122984 -0.90227914
## 7 -1.00401547 1.0511366 1.54783192 -0.90571840 1.55804788
## 8 0.37250149 -2.2951325 0.36548375 -0.98074680 -0.47706564
## 9 0.37784815 -2.1896365 0.57738642 -1.07193157 -0.21050733
## 10 0.47516834 -1.4973383 1.92594459 -1.88157243 1.81617688
## 11 -0.15460428 -1.4874470 0.04500407 -0.56480004 -0.88077329
## 12 -0.75125616 -0.6912640 -0.82286615 1.40572963 0.60861209
## 13 -0.71522346 -0.6123265 -0.68593468 1.22876500 0.95086684
## 14 -0.75371127 -0.5230099 -0.47749143 1.29758114 0.98478126
## 15 -0.68515119 -0.5968216 -0.67456220 1.11594801 1.12060558
## 16 -0.70208367 -0.1179275 -0.65409902 0.13905901 0.22346595
## 17 0.02926667 -1.5816522 -0.82418771 -0.63877930 -0.60848361
## 18 -0.81623681 0.7679641 -0.71900411 -1.66413994 -0.55761297
## 19 -0.83966938 0.4861917 -1.27880250 -1.38952464 -1.31037187
## 20 -0.20134918 -1.4209694 -1.85930318 0.65326621 0.46835677
## 21 -1.32980601 1.1104388 -1.98153818 0.16075938 -0.47401854
## 22 -1.34102037 1.1239204 -1.94642650 0.18949414 -0.49299523
## 23 7.67107700 4.0148949 -1.12616402 -0.80284424 0.80292047 spawndate_pca_data <- spawndate_res_pca$ind$coord[,c(1,2)]
#Creating a data frame containing the earliest spawn date per ECN location, pond, and year, and the co-ordinates of PC1 and PC2
spawndate_glm_data <- cbind(spawndate_pca_data,spawndate_nona)
spawndate_glm_data <- as.data.frame(spawndate_glm_data)
count(spawndate_glm_data) #23 data points## n
## 1 23## Standardise
spawndate_glm_data$z_year<-as.numeric(scale(spawndate_glm_data$year))
spawndate_glm_data$z_tmax<-as.numeric(scale(spawndate_glm_data$tmax))Correlations
###Model
spawndate_glm_data %>% cor_test(.,spawn,Dim.1,method="spearman")## # A tibble: 1 × 6
## var1 var2 cor statistic p method
## <chr> <chr> <dbl> <dbl> <dbl> <chr>
## 1 spawn Dim.1 -0.45 2936. 0.031 Spearman spawndate_glm_data %>% cor_test(.,spawn,Dim.2,method="spearman")## # A tibble: 1 × 6
## var1 var2 cor statistic p method
## <chr> <chr> <dbl> <dbl> <dbl> <chr>
## 1 spawn Dim.2 0.18 1668. 0.422 SpearmanMultivariate Model
#POND ID
with(spawndate_df,table(location,LCODE))## LCODE
## location 1 2 3 4 5 6
## T01 17 0 0 0 0 0
## T02 2 21 0 0 0 0
## T03 18 1 0 0 0 0
## T04 1 7 2 10 18 14
## T05 22 0 0 0 0 0
## T06 13 0 0 0 0 0
## T11 16 13 0 0 0 0 spawndate_df$pond_id<-paste(spawndate_df$location,spawndate_df$LCODE,sep="_")
spawndate_df_subset<-subset(spawndate_df,!is.nan(spawn) & !is.nan(tmax) & !is.na(tmax))
nrow(spawndate_df_subset) #172## [1] 172 table(spawndate_df_subset$location)##
## T01 T02 T03 T04 T05 T06 T11
## 17 23 19 52 19 13 29#library(brms)
bf_spawn<-bf(spawn ~ year + ar(p=1) + (1|p|pond_id)) + gaussian()
bf_spawn_temp<-bf(tmax ~ year + ar(p=1) + (1|p|pond_id)) + gaussian()
#bf_temp_min<-bf(tmin ~ year + (1|p|location)) + gaussian()
###TMAX + TMIN
fit_spawn1 <- brm(
bf_spawn + bf_spawn_temp + set_rescor(TRUE),
data = spawndate_df, chains = 4, cores = 4,
control = list(adapt_delta = 0.99)
)## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
## ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling summary(fit_spawn1)## Family: MV(gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: spawn ~ year + ar(p = 1) + (1 | p | pond_id)
## tmax ~ year + ar(p = 1) + (1 | p | pond_id)
## Data: spawndate_df (Number of observations: 172)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Correlation Structures:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ar_spawn[1] 0.47 0.08 0.32 0.61 1.00 4727 2948
## ar_tmax[1] 0.82 0.05 0.73 0.92 1.00 3258 2431
##
## Multilevel Hyperparameters:
## ~pond_id (Number of levels: 15)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(spawn_Intercept) 23.75 4.74 16.41 34.92 1.00
## sd(tmax_Intercept) 1.70 0.35 1.16 2.55 1.00
## cor(spawn_Intercept,tmax_Intercept) -0.55 0.18 -0.84 -0.12 1.00
## Bulk_ESS Tail_ESS
## sd(spawn_Intercept) 1280 2060
## sd(tmax_Intercept) 1055 2012
## cor(spawn_Intercept,tmax_Intercept) 1277 1641
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## spawn_Intercept 611.37 478.39 -370.24 1508.68 1.00 5209 2897
## tmax_Intercept 9.38 84.72 -155.15 172.16 1.00 3165 2143
## spawn_year -0.27 0.24 -0.72 0.22 1.00 5226 2897
## tmax_year -0.00 0.04 -0.08 0.08 1.00 3158 2144
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_spawn 9.50 0.54 8.52 10.59 1.00 4062 2832
## sigma_tmax 0.64 0.04 0.57 0.72 1.00 4888 3051
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(spawn,tmax) -0.10 0.08 -0.26 0.06 1.00 4535 2875
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1). #Check
pp_check(fit_spawn1,resp='spawn')## Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_check(fit_spawn1,resp='tmax')## Using 10 posterior draws for ppc type 'dens_overlay' by default.
Multivariate Predictions
fit_spawn1_resid<-resid(fit_spawn1)
spawn_resid_plotdata<- data.frame(cbind(fit_spawn1_resid[,,'spawn'][,c(1,3,4)],fit_spawn1_resid[,,'tmax'][,c(1,3,4)]))
colnames(spawn_resid_plotdata)<-c("spawn","spawn2","spawn97","Temp","Temp2","Temp97")
spawn_resid_plotdata$pond_id<-spawndate_df_subset$pond_id
spawn_resid_plotdata$location<-substr(spawn_resid_plotdata$pond_id,1,3)
spawn_plot1<-ggplot(spawn_resid_plotdata,aes(x=Temp,y=spawn))+ theme_bw() + geom_smooth(method="lm")
spawn_plot2<- spawn_plot1 + geom_errorbar(aes(ymin=spawn2,ymax=spawn97,width=0.15),color="gray77") + geom_errorbar(aes(xmin=Temp2,xmax=Temp97,width=2),color="gray77") + geom_point(data=spawn_resid_plotdata,shape=21,size=5,aes(fill=location),color="white") + labs(x="Mean Winter Maximum Temperature Residual",y="Spawn Date Residual",fill="Pond ID") + plotopts + theme(legend.position = "top") + scale_fill_brewer(palette = "Set2") + annotate("text",label="cor = -0.1 [-0.26,0.05]",x=-2,y=55,size=5)
spawn_plot2## `geom_smooth()` using formula = 'y ~ x'
#Mean Shift in Spawn (Days) per 1 deg rise
ggeffects::ggpredict(glm(spawn~Temp,data=spawn_resid_plotdata),terms=list(Temp=1))## # Predicted values of spawn
##
## Temp | Predicted | 95% CI
## ------------------------------
## 1 | -1.58 | -4.18, 1.02Site Level Posterior Plot
spawn_fit1_site<-ranef(fit_spawn1)$pond_id
spawn_site_plotdata<- data.frame(cbind(spawn_fit1_site[,,'spawn_Intercept'][,c(1,3,4)],spawn_fit1_site[,,'tmax_Intercept'][,c(1,3,4)]))
colnames(spawn_site_plotdata)<-c("spawn","spawn2","spawn97","Temp","Temp2","Temp97")
spawn_site_plotdata$location<-substr(rownames(spawn_site_plotdata),1,3)
spawn_site_plot1<-ggplot(spawn_site_plotdata,aes(x=Temp,y=spawn))+ theme_bw() + geom_smooth(method="lm")
spawn_site_plot2<- spawn_site_plot1 + geom_errorbar(aes(ymin=spawn2,ymax=spawn97,width=0.15),color="gray55") + geom_errorbar(aes(xmin=Temp2,xmax=Temp97,width=5),color="gray55") + geom_point(data=spawn_site_plotdata,shape=21,size=5,aes(fill=location),color="white") + labs(x="Mean Maximum Temperature Residual",y="Spawn Date Residual") + plotopts + scale_fill_manual(values=brewer.pal(8,"Set2")[c(1:6,8)]) + annotate("text",label="cor = -0.55 [-0.83,-0.13]",x=-1.5,y=-50,size=6) + theme(legend.position = "top") + labs(fill="")
spawn_site_plot2## `geom_smooth()` using formula = 'y ~ x'
Graphs
Tests
spawndate_long<- spawndate_df %>% dplyr::select(spawn,forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain) %>% pivot_longer(.,cols=c(forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain),names_to="variable",values_to="traitval")
spawndate_corr<-spawndate_long %>% group_by(variable) %>% cor_test(traitval,spawn,method="spearman")
spawndate_corr<-mutate(spawndate_corr,p.adj=p.adjust(p,"BH"))
spawndate_corr$Variable<-label_df$new[match(spawndate_corr$variable,label_df$old)]
spawndate_corr%>% dplyr::select(Variable,cor,statistic,p,p.adj) %>% arrange(p.adj)## # A tibble: 10 × 5
## Variable cor statistic p p.adj
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Grassland Area 0.74 185852. 4.27e-29 4.27e-28
## 2 Arable Area -0.65 1165636. 1.96e-20 9.8 e-20
## 3 Urban Area -0.62 1446777. 6 e-20 2 e-19
## 4 Maximum Temperature -0.61 1365614. 6.24e-19 1.56e-18
## 5 Minimum Temperature -0.59 1350776. 1.06e-17 2.12e-17
## 6 Other Area -0.52 1361490. 9.60e-14 1.6 e-13
## 7 Air Frost 0.53 402149. 1.3 e-13 1.86e-13
## 8 Forest Area -0.53 1083373. 4.62e-13 5.78e-13
## 9 Freshwater Area -0.49 1333948. 3.97e-12 4.41e-12
## 10 Rain 0.13 774558. 7.97e- 2 7.97e- 2Graphs
#Arable area
spawndate_plot1 <- ggplot(spawndate_df, aes(arable_area/1000000, spawn)) +
xlab(bquote("Arable area" (km^2))) + ylab("Spawn Date \n(Days 1st Jan)") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2") + plotopts+ theme(legend.position = "top") + annotate("text",label="cor = -0.65 ***" ,x=7.5,y=100,size=7)
spawndate_plot1## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Grassland
spawndate_plot2 <- ggplot(spawndate_df, aes(grassland_area/1000000, spawn)) +
xlab(bquote("Grassland area" (km^2))) +
ylab("Spawn Date \n(Days 1st Jan)") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2")+ plotopts+ theme(legend.position = "top") + annotate("text",label="cor = 0.74 ***" ,x=8.5,y=100,size=7)
spawndate_plot2 ## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Max Temp
spawndate_plot3 <- ggplot(spawndate_df, aes(tmax, spawn)) +
xlab("Mean Maximum Winter Temperature (\u00B0C)") + ylab("Spawn Date \n(Days 1st Jan)") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = -0.61 ***" ,x=5,y=25,size=7)
spawndate_plot3## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Air Frost
spawndate_plot4 <- ggplot(spawndate_df, aes(af, spawn)) +
xlab("Air Frost Days") + ylab("Spawn Date \n(Days 1st Jan)") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = 0.53 ***" ,x=17.5,y=25,size=7)
spawndate_plot4## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Max Temp Versin 2
spawndate_plot3_v2 <- ggplot(spawndate_df, aes(tmax, spawn,group=location)) +
xlab("Mean Maximum Winter Temperature (\u00B0C)") + ylab("Spawn Date \n(Days 1st Jan)") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = -0.61 ***" ,x=5,y=25,size=7)
spawndate_plot3_v2## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Air Frost Version 2
spawndate_plot4 <- ggplot(spawndate_df, aes(af, spawn)) +
xlab("Air Frost Days") + ylab("Spawn Date \n(Days 1st Jan)") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
theme_bw(base_size = 15) + labs(fill="Pond ID") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = 0.53 ***" ,x=17.5,y=25,size=7)
spawndate_plot4## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
Combined Plot
Spawn All Data
# spawn_multivariateplot<-plot_grid(spawn_siteplot2,spawn_plot2,nrow=1,labels="AUTO",label_size=25)
# ggsave2('Spawn Multivariate Plot.pdf',spawn_multivariateplot,width=35,height=20,units="cm")Individual Predictors
#Save PCA Plot
ggsave2('Spawn PCA PLot.pdf',spawndate_rotationplot,width=10,height=10,units="cm")
#Individual Predictors
spawndate_combinedplot<-plot_grid(spawndate_plot1,spawndate_plot2,spawndate_plot3,spawndate_plot4,labels=c("C","D","E","F"),nrow=2,label_size = 25)## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).#Stich On Site Specific Temp Trends
spawndate_combinedplot2<-plot_grid(spawn_siteplot2,spawn_site_plot2,ncol=2,labels = c("A","B"),label_size = 25)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'#Stich On Bayesian Model
spawndate_combinedplot3<-plot_grid(spawndate_combinedplot2,spawndate_combinedplot,nrow=2,rel_heights = c(1.5,2))
#Save
ggsave2('Fig_2_Spawn Grid Plot.pdf',spawndate_combinedplot3,width=30,height=45,units="cm")HATCHDATE
Mean Trend Over Time
# Correlations
hatchdate_df %>% group_by(location) %>% cor_test(hatch,year,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))## # A tibble: 7 × 10
## location var1 var2 cor statistic p conf.low conf.high method p.adj
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 T01 hatch year -0.23 -0.535 0.616 -0.839 0.631 Pears… 0.677
## 2 T02 hatch year -0.36 -1.75 0.0944 -0.671 0.0646 Pears… 0.165
## 3 T03 hatch year 0.57 2.30 0.0421 0.0274 0.853 Pears… 0.0982
## 4 T04 hatch year -0.4 -2.87 0.00642 -0.621 -0.121 Pears… 0.0449
## 5 T05 hatch year 0.097 0.423 0.677 -0.350 0.507 Pears… 0.677
## 6 T06 hatch year 0.45 1.41 0.197 -0.256 0.840 Pears… 0.276
## 7 T11 hatch year -0.44 -2.41 0.0239 -0.708 -0.0656 Pears… 0.0836#Models
hatch_mods<-hatchdate_df %>% group_by(location) %>% do(tidy(lm(hatch~year,data=.))) %>% filter(term=="year")
hatch_mods$p.adj<-p.adjust(p=hatch_mods$p.value,"BH")
hatch_mods## # A tibble: 7 × 7
## # Groups: location [7]
## location term estimate std.error statistic p.value p.adj
## <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 T01 year -0.452 0.845 -0.535 0.616 0.677
## 2 T02 year -0.576 0.329 -1.75 0.0944 0.165
## 3 T03 year 1.32 0.574 2.30 0.0421 0.0982
## 4 T04 year -0.764 0.267 -2.87 0.00642 0.0449
## 5 T05 year 0.147 0.346 0.423 0.677 0.677
## 6 T06 year 0.752 0.534 1.41 0.197 0.276
## 7 T11 year -1.01 0.419 -2.41 0.0239 0.0835 #
#Sample Size
hatchdate_df %>% filter(!is.na(hatch)) %>% nrow()## [1] 145hatchmeans<- hatchdate_df %>% group_by(year,location) %>% reframe(hatch=mean(hatch,na.rm=T))
hatch_site_plot1<-ggplot(hatchdate_df,aes(x=year,y=hatch))+ geom_smooth(method="lm") + geom_point(shape=21,size=5,aes(fill=location)) + theme_bw() + facet_wrap(.~location) + labs(x="Year",y="Hatch Date (Days 1st Jan)") + plotopts + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2") + guides(fill="none")
hatch_site_plot1## `geom_smooth()` using formula = 'y ~ x'
#Annotations
hatch_ann_text <- data.frame(year = 2005,hatch = 60,lab = "cor = -0.4 \np.adj = 0.045",
location = "T04")
hatch_site_plot2<-hatch_site_plot1 + geom_text(data = hatch_ann_text,label = hatch_ann_text$lab,size=6)
hatch_site_plot2## `geom_smooth()` using formula = 'y ~ x'
Mean Trend with Temperature
hatch_tmax_plot1<-ggplot(hatchdate_df)+ geom_smooth(method="lm",aes(x=tmax,y=hatch,group=location)) + geom_point(shape=21,size=5,aes(fill=location,x=tmax,y=hatch)) + guides(fill="none") + labs(x="Mean Spawning Temperature",y="Hatch Date (Days 1st Jan)") + plotopts + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2")
hatch_tmax_plot1## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Correlations
spawndate_df %>% group_by(location) %>% cor_test(spawn,tmax,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))## # A tibble: 7 × 10
## location var1 var2 cor statistic p conf.low conf.high method p.adj
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 T01 spawn tmax -0.31 -1.28 0.219 -0.691 0.196 Pearson 0.307
## 2 T02 spawn tmax -0.45 -2.33 0.0298 -0.729 -0.0508 Pearson 0.0695
## 3 T03 spawn tmax -0.19 -0.804 0.433 -0.594 0.288 Pearson 0.449
## 4 T04 spawn tmax -0.34 -2.56 0.0135 -0.561 -0.0746 Pearson 0.0695
## 5 T05 spawn tmax -0.52 -2.52 0.0218 -0.789 -0.0892 Pearson 0.0695
## 6 T06 spawn tmax -0.37 -1.31 0.217 -0.764 0.230 Pearson 0.307
## 7 T11 spawn tmax -0.15 -0.768 0.449 -0.487 0.233 Pearson 0.449PCA
#Removing sun
hatchdate_nosun <- hatchdate_df %>% dplyr::select(-sun)
#Removing all rows containing NA values
hatchdate_nona <- na.omit(hatchdate_nosun)
#Removing response variable
hatchdate_noresp <- hatchdate_nona %>% dplyr::select(-hatch)
#Removing location variables
hatchdate_for_pca <- hatchdate_noresp %>% dplyr::select(-location)
hatchdate_for_pca <- hatchdate_for_pca %>% dplyr::select(-LCODE)
#Variable Rename
hatchdatepcalabs<-label_df$new[match(colnames(hatchdate_for_pca),label_df$old)]
#hatchdatepcalabs[c(2,3)]<-c("Pre-Spawn Ammonium","Pre-Spawn Nitrate")
colnames(hatchdate_for_pca)<-hatchdatepcalabs
#Performing Principal Components Analysis
hatchdate_res_pca <- PCA( hatchdate_for_pca, graph = FALSE)
#Contributions
dimdesc(hatchdate_res_pca, axes = 1:2)## $Dim.1
##
## Link between the variable and the continuous variables (R-square)
## =================================================================================
## correlation p.value
## Arable Area 0.9525272 8.620719e-12
## Urban Area 0.9256264 6.858454e-10
## Maximum Temperature 0.8416292 9.158849e-07
## Minimum Temperature 0.7354440 9.621395e-05
## Spawn Nitrate 0.6012114 3.082123e-03
## Air Frost -0.6891667 3.889722e-04
## Grassland Area -0.9254699 6.999493e-10
##
## $Dim.2
##
## Link between the variable and the continuous variables (R-square)
## =================================================================================
## correlation p.value
## Other Area 0.8670282 1.781389e-07
## Freshwater Area 0.8210231 2.842221e-06
## Minimum Temperature 0.6394147 1.355243e-03
## Forest Area -0.5282852 1.149194e-02
## Air Frost -0.6774766 5.325217e-04#Producing a rotation plot
hatch_rotationplot<-fviz_pca_var(hatchdate_res_pca, col.var = "black",
repel = TRUE
) +
labs(x = "PC1", y = "PC2", title = "") + guides(color="none")Data
#Creating a data frame containing the co-ordinates of PC1 and PC2
hatchdate_res_pca$ind$coord## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 2.1624428 0.25964591 -1.2696770 0.38552765 -0.838354052
## 2 -1.4517261 -0.51039504 0.5089993 -2.09556132 -0.670677701
## 3 2.0207646 -0.45675381 -0.7301879 -0.78668106 0.756124000
## 4 -0.2503202 0.68325087 -1.8213885 -2.02786186 0.437144289
## 5 -0.7179078 0.44545375 -0.9770840 -0.97576641 -1.166821965
## 6 3.0765488 0.30719018 -2.5877631 -1.05572833 2.394646993
## 7 -0.4474291 3.50602841 0.6953934 0.37425726 0.030461439
## 8 -1.0329187 0.06097587 -1.6724709 -0.36589037 -1.573565459
## 9 -1.0421476 0.05873357 -1.6125910 -0.31524971 -1.665392671
## 10 -1.4359662 -1.56826957 -0.4754190 0.06326624 0.888455964
## 11 6.1504394 -2.06732928 3.5165970 -1.48393673 -0.699559162
## 12 -0.3945104 2.78829342 1.1483096 -0.06023258 0.959020762
## 13 -0.4389931 2.77865336 1.4257118 0.17763260 0.531418226
## 14 -0.8021923 2.30953440 1.0896670 0.84135802 -0.108996229
## 15 -0.7959737 2.30895969 1.0693540 0.81834754 -0.073916461
## 16 -1.7241263 -2.07002364 0.1277186 1.08087596 0.393895395
## 17 -2.7973736 -2.98795026 0.8342562 0.73967544 1.212208737
## 18 -2.3445280 -1.50136397 1.0276819 -0.85785737 0.187863205
## 19 -2.3948673 -1.51089401 1.3283574 -0.59602348 -0.278317509
## 20 -1.3346020 -1.52371074 -0.5712024 1.89469024 -0.006008684
## 21 3.2451164 -0.33102203 -1.1280507 2.34477954 -0.606378753
## 22 2.7502704 -0.97900708 0.0737884 1.90037875 -0.103250362 hatchdate_pca_data <- hatchdate_res_pca$ind$coord[,c(1,2)]
#Creating a data frame containing the earliest hatch date per ECN location, pond, and year, and the co-ordinates of PC1 and PC2
hatchdate_glm_data <- cbind(hatchdate_pca_data,hatchdate_nona)
colnames(hatchdate_glm_data) [6]<- c("hatchdate")
hatchdate_glm_data <- as.data.frame(hatchdate_glm_data)
count(hatchdate_glm_data) #25 data points## n
## 1 22 table(hatchdate_glm_data$location)##
## T01 T02 T03 T04 T05 T06 T11
## 1 4 0 7 5 0 5 nrow(hatchdate_glm_data)## [1] 22Correlations
###Model
hatchdate_glm_data %>% cor_test(.,hatchdate,Dim.1,method="spearman")## # A tibble: 1 × 6
## var1 var2 cor statistic p method
## <chr> <chr> <dbl> <dbl> <dbl> <chr>
## 1 hatchdate Dim.1 -0.7 3009. 0.000293 Spearman hatchdate_glm_data %>% cor_test(.,hatchdate,Dim.2,method="spearman")## # A tibble: 1 × 6
## var1 var2 cor statistic p method
## <chr> <chr> <dbl> <dbl> <dbl> <chr>
## 1 hatchdate Dim.2 0.033 1712. 0.883 SpearmanMultivariate Model
#POND ID
with(hatchdate_df,table(location,LCODE))## LCODE
## location 1 2 4 5 6 7
## T01 7 0 0 0 0 0
## T02 2 21 0 0 0 0
## T03 12 1 0 0 0 0
## T04 0 5 8 18 13 1
## T05 21 0 0 0 0 0
## T06 10 0 0 0 0 0
## T11 15 11 0 0 0 0 hatchdate_df$pond_id<-paste(hatchdate_df$location,hatchdate_df$LCODE,sep="_")
hatchdate_df_subset<-subset(hatchdate_df,!is.nan(hatch) & !is.nan(tmax) & !is.na(tmax))
nrow(hatchdate_df_subset) #143## [1] 143#library(brms)
bf_hatch<-bf(hatch ~ year + ar(p=1) + (1|p|pond_id)) + gaussian()
bf_hatch_temp<-bf(tmax ~ year + ar(p=1) + (1|p|pond_id)) + gaussian()
#bf_temp_min<-bf(tmin ~ year + (1|p|location)) + gaussian()
###TMAX + TMIN
fit_hatch1 <- brm(
bf_hatch + bf_hatch_temp + set_rescor(TRUE),
data = hatchdate_df_subset, chains = 4, cores = 4,
control = list(adapt_delta = 0.99)
)## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
## ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling summary(fit_hatch1)## Family: MV(gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: hatch ~ year + ar(p = 1) + (1 | p | pond_id)
## tmax ~ year + ar(p = 1) + (1 | p | pond_id)
## Data: hatchdate_df_subset (Number of observations: 143)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Correlation Structures:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ar_hatch[1] 0.34 0.08 0.18 0.50 1.00 3959 2917
## ar_tmax[1] 0.66 0.07 0.52 0.80 1.00 3854 2810
##
## Multilevel Hyperparameters:
## ~pond_id (Number of levels: 14)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(hatch_Intercept) 22.26 4.55 15.28 32.95 1.00
## sd(tmax_Intercept) 1.49 0.33 1.00 2.27 1.00
## cor(hatch_Intercept,tmax_Intercept) -0.47 0.21 -0.80 0.02 1.00
## Bulk_ESS Tail_ESS
## sd(hatch_Intercept) 1421 2085
## sd(tmax_Intercept) 1286 1864
## cor(hatch_Intercept,tmax_Intercept) 1437 1851
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## hatch_Intercept 758.96 437.35 -111.00 1618.20 1.00 4250 2971
## tmax_Intercept -44.80 50.75 -144.45 49.17 1.00 3606 2464
## hatch_year -0.33 0.22 -0.76 0.11 1.00 4256 2989
## tmax_year 0.03 0.03 -0.02 0.08 1.00 3609 2464
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_hatch 9.82 0.62 8.71 11.12 1.00 3957 2926
## sigma_tmax 0.60 0.04 0.53 0.68 1.00 4258 2895
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(hatch,tmax) -0.37 0.08 -0.52 -0.20 1.00 3653 2898
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1). bayes_R2(fit_hatch1)## Estimate Est.Error Q2.5 Q97.5
## R2hatch 0.8409414 0.01221426 0.8126497 0.8600146
## R2tmax 0.8497104 0.01190528 0.8226606 0.8682156 #Check
# pp_check(fit_spawn1,resp='highestpercdeadlogit')
# pp_check(fit_spawn1,resp='tmax')Multivariate Predictions
fit_hatch1_resid<-resid(fit_hatch1)
hatch_resid_plotdata<- data.frame(cbind(fit_hatch1_resid[,,'hatch'][,c(1,3,4)],fit_hatch1_resid[,,'tmax'][,c(1,3,4)]))
colnames(hatch_resid_plotdata)<-c("hatch","hatch2","hatch97","Temp","Temp2","Temp97")
hatch_resid_plotdata$pond_id<-hatchdate_df_subset$pond_id
hatch_resid_plotdata$location<-hatchdate_df_subset$location
hatch_plot1<-ggplot(hatch_resid_plotdata,aes(x=Temp,y=hatch))+ theme_bw() + geom_smooth(method="lm")
hatch_plot2<- hatch_plot1 + geom_errorbar(aes(ymin=hatch2,ymax=hatch97,width=0.15),color="gray77") + geom_errorbar(aes(xmin=Temp2,xmax=Temp97,width=2),color="gray77") + geom_point(data=hatch_resid_plotdata,shape=21,size=5,aes(fill=location),color="white") + labs(x="Mean Maximum Spawning \nTemperature Residual",y="Hatch Date Residual",fill="Location") + plotopts + theme(legend.position = "top") + scale_fill_brewer(palette="Set2") + annotate("text",x=-1.8,y=-25,label="cor = -0.37 \n[-0.52, -0.22]",size=6)
hatch_plot2## `geom_smooth()` using formula = 'y ~ x'
#Shift in Hatching (Days) Per 1 deg rise in temperature
ggeffects::ggpredict(glm(hatch~Temp,data=hatch_resid_plotdata),terms=list(Temp=1))## # Predicted values of hatch
##
## Temp | Predicted | 95% CI
## -------------------------------
## 1 | -6.24 | -9.10, -3.38Graphs
Tests
hatchdate_long<- hatchdate_df %>% dplyr::select(hatch,forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain) %>% pivot_longer(.,cols=c(forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain),names_to="variable",values_to="traitval")
hatchdate_corr<-hatchdate_long %>% group_by(variable) %>% cor_test(traitval,hatch,method="spearman")
hatchdate_corr<-mutate(hatchdate_corr,p.adj=p.adjust(p,"BH"))
hatchdate_corr$Variable<-label_df$new[match(hatchdate_corr$variable,label_df$old)]
hatchdate_corr%>% dplyr::select(Variable,cor,statistic,p,p.adj) %>% arrange(p.adj)## # A tibble: 10 × 5
## Variable cor statistic p p.adj
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Grassland Area 0.74 106183. 9.12e-25 9.12e-24
## 2 Urban Area -0.69 857599. 1.20e-21 6 e-21
## 3 Arable Area -0.68 687733. 1.91e-19 6.37e-19
## 4 Minimum Temperature -0.61 784169. 6.98e-16 1.74e-15
## 5 Air Frost 0.54 224812. 3.92e-12 7.84e-12
## 6 Maximum Temperature -0.51 737671. 5.36e-11 8.93e-11
## 7 Forest Area -0.52 624881. 6.97e-11 9.96e-11
## 8 Other Area -0.4 713472. 4.59e- 7 5.74e- 7
## 9 Freshwater Area -0.4 711331. 6.19e- 7 6.88e- 7
## 10 Rain 0.18 398811. 2.99e- 2 2.99e- 2Graphs
#Arable area
hatchdate_plot1 <- ggplot(hatchdate_df, aes(arable_area/1000000, hatch)) +
xlab(bquote("Arable area" (km^2))) + ylab("Hatch Date (Days 1st Jan)") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + plotopts+ theme(legend.position = "top") + annotate("text",label="cor = -0.68 ***" ,x=7.5,y=125,size=7)
hatchdate_plot1## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
#T MAX
hatchdate_plot2 <- ggplot(hatchdate_df, aes(tmax, hatch)) +
xlab("Mean Maximum Spawning Temperature") +
ylab("Hatch Date (Days 1st Jan)") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2")+ plotopts+ theme(legend.position = "top") + annotate("text",label="cor = -0.51 ***" ,x=8,y=50,size=7)
hatchdate_plot2 ## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Grassland area
hatchdate_plot3 <- ggplot(hatchdate_df, aes(grassland_area/1000000, hatch)) +
xlab(bquote("Grassland area" (km^2))) + ylab("Hatch Date (Days 1st Jan)") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = 0.74 ***" ,x=9,y=120,size=7)
hatchdate_plot3## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Air Frost
hatchdate_plot4 <- ggplot(hatchdate_df, aes(af, hatch)) +
xlab("Air Frost Days per Month") +
ylab("Hatch date (Days 1st Jan)") +
geom_smooth(method = "lm",formula=y~poly(x,3)) + geom_point(shape=21,size=5,aes(fill=location),color="white") +
theme_bw(base_size = 15) + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top") + annotate("text",label="cor = 0.54 ***" ,x=10,y=50,size=7)
hatchdate_plot4## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Urban area
hatchdate_plot5 <- ggplot(hatchdate_df, aes(urban_area/1000000, hatch)) +
xlab(bquote("Urban area" (km^2))) + ylab("Hatch Date (Days 1st Jan)") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location),color="white") +
theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + plotopts + theme(legend.position = "top")
hatchdate_plot5## `geom_smooth()` using formula = 'y ~ x'
Individual Predictors
#Save PCA Plot
ggsave2('Hatch PCA PLot.pdf',hatch_rotationplot,width=10,height=10,units="cm")
#Individual Predictors
hatch_combinedplot<-plot_grid(hatchdate_plot1,hatchdate_plot3,hatchdate_plot2,hatchdate_plot4,labels=c("C","D","E","F"),nrow=2,label_size = 25)## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).#Stich On Site Specific Temp Trends
hatchdate_combinedplot2<-plot_grid(hatch_site_plot2,hatch_plot2,ncol=2,labels = c("A","B"),label_size = 25)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'#Stich On Bayesian Model
hatchdate_combinedplot3<-plot_grid(hatchdate_combinedplot2,hatch_combinedplot,nrow=2,rel_heights = c(1.5,2))
#Save
ggsave2('Fig_3_Hatch Grid Plot.pdf',hatchdate_combinedplot3,width=30,height=45,units="cm")PERCENT DEAD
Mean Trend Over Time
percdead_df %>% group_by(year,location) %>% reframe(percdead=mean(percdead,na.rm=T)) ## # A tibble: 119 × 3
## year location percdead
## <int> <fct> <dbl>
## 1 1994 T01 32.9
## 2 1994 T02 15
## 3 1994 T03 0.167
## 4 1994 T04 4.40
## 5 1994 T05 2.14
## 6 1994 T06 0
## 7 1995 T01 33.9
## 8 1995 T02 2
## 9 1995 T03 2.08
## 10 1995 T04 17.0
## # ℹ 109 more rowspercdead_df %>% group_by(location) %>% reframe(min=min(year,na.rm=T),max=max(year,na.rm=T)) ## # A tibble: 7 × 3
## location min max
## <fct> <int> <int>
## 1 T01 1994 2010
## 2 T02 1994 2014
## 3 T03 1994 2006
## 4 T04 1994 2015
## 5 T05 1994 2014
## 6 T06 1994 2010
## 7 T11 1995 2011hatchdate_df %>% group_by(location) %>% reframe(min=min(year,na.rm=T),max=max(year,na.rm=T)) ## # A tibble: 7 × 3
## location min max
## <fct> <int> <int>
## 1 T01 1995 2005
## 2 T02 1994 2014
## 3 T03 1994 2006
## 4 T04 1994 2015
## 5 T05 1994 2015
## 6 T06 1994 2010
## 7 T11 1996 2012spawndate_df %>% group_by(location) %>% reframe(min=min(year,na.rm=T),max=max(year,na.rm=T)) ## # A tibble: 7 × 3
## location min max
## <fct> <int> <int>
## 1 T01 1994 2010
## 2 T02 1994 2014
## 3 T03 1994 2011
## 4 T04 1994 2015
## 5 T05 1994 2015
## 6 T06 1994 2010
## 7 T11 1996 2012#Subset
percdead_df_subset<-subset(percdead_df,percdead>=0)
#SampleSize
percdead_df %>% filter(!is.nan(percdead)) %>% nrow()## [1] 178# Correlations
percdead_df %>% filter(!is.nan(percdead)) %>% group_by(location) %>% cor_test(percdead,year,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))## # A tibble: 7 × 10
## location var1 var2 cor statistic p conf.low conf.high method p.adj
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 T01 perc… year -0.043 -0.161 0.875 -5.27e-1 0.463 Pears… 0.882
## 2 T02 perc… year 0.42 2.08 0.0505 2.46e-4 0.716 Pears… 0.177
## 3 T03 perc… year -0.34 -1.23 0.241 -7.35e-1 0.238 Pears… 0.506
## 4 T04 perc… year 0.018 0.149 0.882 -2.23e-1 0.257 Pears… 0.882
## 5 T05 perc… year 0.26 1.10 0.289 -2.31e-1 0.651 Pears… 0.506
## 6 T06 perc… year -0.17 -0.556 0.59 -6.56e-1 0.424 Pears… 0.826
## 7 T11 perc… year 0.54 3.27 0.00301 2.09e-1 0.760 Pears… 0.0211 # percdead_df_subset %>% filter(!is.nan(percdead)) %>% group_by(pond_id) %>% cor_test(percdead,year,method="spearman") %>% mutate(p.adj=p.adjust(p,method="BH"))
#Plot By Location
percdeaad_siteplot1<-ggplot(percdead_df,aes(x=year,y=percdead))+ geom_smooth(method="lm") + geom_point(shape=21,size=5,aes(fill=location)) + theme_bw() + facet_wrap(.~location) + guides(fill="none") + plotopts + labs(x="Year",y="Percetage Spawn Dead") + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2")
percdeaad_siteplot1## `geom_smooth()` using formula = 'y ~ x'
#Plot By Pond ID
# percdeaad_siteplot2<-ggplot(percdead_df_subset,aes(x=year,y=percdead))+ geom_smooth(method="lm") + geom_point(shape=21,size=5,aes(fill=location)) + theme_bw() + facet_wrap(.~pond_id) + guides(fill="none") + plotopts + labs(x="Year",y="Percetage Spawn Dead") + theme(axis.text.x = element_text(angle=45,hjust=1)) + scale_fill_brewer(palette = "Set2")
# percdeaad_siteplot2
#Annotations
ann_text_percdead <- data.frame(year = 2001,percdead = 75,lab = "cor = 0.54 \np.adj = 0.02",
location = "T11")
percdeaad_siteplot1_ann<-percdeaad_siteplot1 + geom_text(data = ann_text_percdead,label = ann_text_percdead$lab,size=5)
percdeaad_siteplot1_ann ## `geom_smooth()` using formula = 'y ~ x'
Mean Trend with Temperature
dead_tmax_plot1<-ggplot(percdead_df)+ geom_smooth(method="lm",aes(x=af,y=percdead,group=location)) + geom_point(shape=21,size=5,aes(fill=location,x=af,y=percdead)) + guides(fill="none") + labs(x="Mean Spawning Temperature",y="Percent Dead Spawn") + plotopts + scale_fill_brewer(palette = "Set2")
dead_tmax_plot1## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Correlations
spawndate_df %>% group_by(location) %>% cor_test(spawn,tmax,method="pearson") %>% mutate(p.adj=p.adjust(p,method="BH"))## # A tibble: 7 × 10
## location var1 var2 cor statistic p conf.low conf.high method p.adj
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 T01 spawn tmax -0.31 -1.28 0.219 -0.691 0.196 Pearson 0.307
## 2 T02 spawn tmax -0.45 -2.33 0.0298 -0.729 -0.0508 Pearson 0.0695
## 3 T03 spawn tmax -0.19 -0.804 0.433 -0.594 0.288 Pearson 0.449
## 4 T04 spawn tmax -0.34 -2.56 0.0135 -0.561 -0.0746 Pearson 0.0695
## 5 T05 spawn tmax -0.52 -2.52 0.0218 -0.789 -0.0892 Pearson 0.0695
## 6 T06 spawn tmax -0.37 -1.31 0.217 -0.764 0.230 Pearson 0.307
## 7 T11 spawn tmax -0.15 -0.768 0.449 -0.487 0.233 Pearson 0.449Multivariate model
#POND ID
#with(percdead_df,table(location,LCODE))
percdead_df$pond_id<-paste(percdead_df$location,percdead_df$LCODE,sep="_")
nrow(percdead_df_subset) ### [1] 178 pdf_subtable<-table(percdead_df$pond_id)
percdead_df_subset<-subset(percdead_df,pond_id %in% names(pdf_subtable)[pdf_subtable>9])
nrow(percdead_df_subset)## [1] 165#Convert To Logits
library(arm)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:rstatix':
##
## select## The following object is masked from 'package:dplyr':
##
## select## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loading required package: lme4##
## Attaching package: 'lme4'## The following object is masked from 'package:brms':
##
## ngrps##
## arm (Version 1.14-4, built: 2024-4-1)## Working directory is /Users/xavharrison/Library/CloudStorage/Dropbox/Mac/Documents/1. Projects/14. Frog Phenology Kat Oliver/Frog Phenology Statistics##
## Attaching package: 'arm'## The following object is masked from 'package:corrplot':
##
## corrplot percdead_df_subset$percdead[which(percdead_df_subset$percdead==100)]<-99.99999
percdead_df_subset$percdead[which(percdead_df_subset$percdead==0)]<-0.000001
percdead_df_subset$percdead_logit<-logit(percdead_df_subset$percdead/100)
## ALL DATA, TMAX
bf_dead<-bf(percdead_logit ~ year + ar(p=1) + (1|p|pond_id)) + gaussian()
#bf_dead_perc<-bf(percdead ~ year + ar(p=1) + (1|p|pond_id)) + gaussian()
#bf_af<-bf(af ~ 1 + ar(p=1) + (0+year|p|pond_id)) + gaussian()
#bf_temp_min<-bf(tmin ~ 1 + ar(p=1) + (1|p|pond_id)) + gaussian()
bf_temp_max<-bf(tmax ~ year + ar(p=1) + (1|p|pond_id)) + gaussian()
# ### TMIN
# percedead_fit1 <- brm(
# bf_dead_perc + bf_temp_max + set_rescor(TRUE),
# data = percdead_df_subset, chains = 4, cores = 4,
# control = list(adapt_delta = 0.99)
# )
#
# summary(percedead_fit1)
#
# #Check
# pp_check(percedead_fit1,resp='percdead')
# pp_check(percedead_fit1,resp='tmax')
### TMAX
percedead_fit2 <- brm(
bf_dead + bf_temp_max + set_rescor(TRUE),
data = percdead_df_subset, chains = 4, cores = 4,
control = list(adapt_delta = 0.99)
)## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
## ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling summary(percedead_fit2)## Family: MV(gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: percdead_logit ~ year + ar(p = 1) + (1 | p | pond_id)
## tmax ~ year + ar(p = 1) + (1 | p | pond_id)
## Data: percdead_df_subset (Number of observations: 163)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Correlation Structures:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ar_percdeadlogit[1] 0.02 0.09 -0.15 0.19 1.00 5237 2635
## ar_tmax[1] 0.82 0.06 0.71 0.93 1.00 3348 2411
##
## Multilevel Hyperparameters:
## ~pond_id (Number of levels: 11)
## Estimate Est.Error l-95% CI
## sd(percdeadlogit_Intercept) 5.32 1.23 3.50
## sd(tmax_Intercept) 1.32 0.29 0.88
## cor(percdeadlogit_Intercept,tmax_Intercept) -0.80 0.13 -0.96
## u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(percdeadlogit_Intercept) 8.28 1.00 1457 1877
## sd(tmax_Intercept) 2.00 1.00 1251 1869
## cor(percdeadlogit_Intercept,tmax_Intercept) -0.45 1.00 1462 2207
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## percdeadlogit_Intercept -79.58 146.65 -369.37 202.82 1.00 4398
## tmax_Intercept -60.87 75.55 -225.05 74.43 1.00 3525
## percdeadlogit_year 0.04 0.07 -0.10 0.18 1.00 4382
## tmax_year 0.04 0.04 -0.03 0.12 1.00 3519
## Tail_ESS
## percdeadlogit_Intercept 3008
## tmax_Intercept 1913
## percdeadlogit_year 2985
## tmax_year 1913
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_percdeadlogit 4.61 0.26 4.13 5.15 1.00 5159 2761
## sigma_tmax 0.54 0.03 0.48 0.61 1.00 3878 2679
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## rescor(percdeadlogit,tmax) 0.04 0.09 -0.13 0.21 1.00 5128
## Tail_ESS
## rescor(percdeadlogit,tmax) 3130
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Site Level Posterior Plot
percdead_fit2_site<-ranef(percedead_fit2)$pond_id
percdead_site_plotdata<- data.frame(cbind(percdead_fit2_site[,,'percdeadlogit_Intercept'][,c(1,3,4)],percdead_fit2_site[,,'tmax_Intercept'][,c(1,3,4)]))
colnames(percdead_site_plotdata)<-c("percdead","percdead2","percdead97","Temp","Temp2","Temp97")
percdead_site_plotdata$location<-substr(rownames(percdead_site_plotdata),1,3)
percdead_site_plot1<-ggplot(percdead_site_plotdata,aes(x=Temp,y=percdead))+ theme_bw() + geom_smooth(method="lm")
percdead_site_plot2<- percdead_site_plot1 + geom_errorbar(aes(ymin=percdead2,ymax=percdead97,width=0.15),color="gray55") + geom_errorbar(aes(xmin=Temp2,xmax=Temp97,width=0.5),color="gray55") + geom_point(data=percdead_site_plotdata,shape=21,size=5,aes(fill=location),color="white") + labs(x="Mean Maximum Temperature Residual",y="Percentage Dead Spawn Residual",fill="Location") + plotopts + scale_fill_brewer(palette="Set2") + annotate(geom="text",x=-1,y=-10,label="cor = -0.8 \n [-0.96,-0.45]",size=6)
percdead_site_plot2## `geom_smooth()` using formula = 'y ~ x'
PCA
#Removing sun
percdead_nosun <- percdead_df %>% dplyr::select(-sun)
#Removing all rows containing NA values
percdead_nona <- na.omit(percdead_nosun)
nrow(percdead_nona)## [1] 29#Removing response variable
percdead_noresp <- percdead_nona %>% dplyr::select(-percdead)
#Removing location variables
percdead_for_pca <- percdead_noresp %>% dplyr::select(-location,-LCODE,-pond_id)
#percdead_for_pca <- percdead_for_pca %>% dplyr::select(-LCODE)
#head(percdead_for_pca)
#Rename Variables
percdeadpcalabs<-label_df$new[match(colnames(percdead_for_pca),label_df$old)]
colnames(percdead_for_pca)<-percdeadpcalabs
#Performing Principal Components Analysis
percdead_res_pca <- PCA(percdead_for_pca, graph = FALSE)
#Contributions
dimdesc(percdead_res_pca, axes = 1:2)## $Dim.1
##
## Link between the variable and the continuous variables (R-square)
## =================================================================================
## correlation p.value
## Arable Area 0.9778107 7.290119e-20
## Urban Area 0.9655034 2.619852e-17
## Spawn Nitrate 0.8648117 1.450928e-09
## Maximum Temperature 0.8121883 8.833833e-08
## Minimum Temperature 0.5908527 7.387844e-04
## Rain -0.4438891 1.585945e-02
## Air Frost -0.5182267 3.979574e-03
## Grassland Area -0.9636134 5.323038e-17
##
## $Dim.2
##
## Link between the variable and the continuous variables (R-square)
## =================================================================================
## correlation p.value
## Other Area 0.8275528 3.075322e-08
## Freshwater Area 0.7758503 7.625761e-07
## Minimum Temperature 0.7536917 2.358379e-06
## Year -0.3765206 4.409220e-02
## Forest Area -0.3786648 4.280448e-02
## Air Frost -0.7995917 1.958766e-07#Producing a rotation plot
percdead_rotationplot<-fviz_pca_var(percdead_res_pca, col.var = "black",
repel = TRUE
) +
labs(x = "PC1", y = "PC2", title = "") + guides(colour="none") + plotoptsData Processing
#Creating a data frame containing the co-ordinates of PC1 and PC2
percdead_res_pca$ind$coord## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 0.7512994 1.02657507 -1.51268259 1.12409659 -1.79597090
## 2 -1.8179391 -0.62499945 -0.46850191 -1.89489783 -0.92005141
## 3 -1.6494503 -0.67708626 -0.20695827 -1.93082601 -1.19829574
## 4 0.6958885 0.18615190 -1.04519424 0.31550973 -0.78384539
## 5 -0.8130605 0.71999639 -1.78429553 -1.01966290 -0.26910012
## 6 -1.1505938 0.45776866 -1.14647909 -0.67115920 -1.12834228
## 7 1.4341075 1.10269584 -2.22329693 0.90449787 0.45993019
## 8 -1.3637829 3.92787997 0.99105864 -0.01834216 -0.01495594
## 9 3.0312263 -0.59885107 -1.29333748 -0.94173832 1.73855861
## 10 -1.5782410 -0.04932380 -1.51337000 -0.06375137 -0.62033092
## 11 -1.5306969 -0.05647557 -3.56977066 -0.71008298 2.53710954
## 12 -1.5786135 -0.04935421 -1.47285467 -0.05139918 -0.68214031
## 13 4.1095036 -0.44927316 0.43559229 -0.61179465 -0.05574561
## 14 -1.8626460 -1.75986632 -0.04224800 0.83052575 -0.19762166
## 15 4.0829075 -0.49860264 1.13086800 -0.78057490 -0.40147144
## 16 -1.0550380 3.09177794 1.50415174 -0.27695781 0.62702660
## 17 -1.0579306 3.09198789 1.69282376 -0.21865238 0.33837409
## 18 3.3708552 -1.02942816 0.97701422 -0.53933873 -0.08047223
## 19 3.9976351 0.01687874 -0.05966336 0.41406738 -0.01262993
## 20 -1.4353333 2.54179092 1.52547898 0.31284058 0.45205583
## 21 -1.4331228 2.54117462 1.50998641 0.30671756 0.47715233
## 22 3.4496714 -0.88494376 1.15779221 -0.13901098 -0.05575183
## 23 3.0210653 -1.22053408 1.13354547 -0.06159835 0.21793728
## 24 -2.0215268 -2.34070991 0.72577747 1.49505406 0.02999764
## 25 -2.9056193 -3.50928310 1.48500017 0.95991272 0.58068448
## 26 -2.2169766 -2.02808633 1.17699888 -1.13560846 0.51422729
## 27 -2.2215458 -2.02742762 1.38270306 -1.07108150 0.19851673
## 28 -1.8499598 -1.56232549 0.23147536 2.44238246 0.12439102
## 29 1.5979173 0.66189299 -0.72161393 3.03087302 -0.07923592 percdead_pca_data <- percdead_res_pca$ind$coord[,c(1,2)]
#Creating a data frame containing the highest percentage of spawn dead per ECN location, pond, and year, and the co-ordinates of PC1 and PC2
percdead_glm_data <- cbind(percdead_pca_data,percdead_nona)
#colnames(percdead_glm_data)[1:2] <- c("Dim.1", "Dim.2")
percdead_glm_data <- as.data.frame(percdead_glm_data)
count(percdead_glm_data) #33 data points## n
## 1 29Correlations
###Model
percdead_glm_data %>% cor_test(.,percdead,Dim.1,method="spearman")## # A tibble: 1 × 6
## var1 var2 cor statistic p method
## <chr> <chr> <dbl> <dbl> <dbl> <chr>
## 1 percdead Dim.1 -0.49 6065. 0.00648 Spearman percdead_glm_data %>% cor_test(.,percdead,Dim.2,method="spearman")## # A tibble: 1 × 6
## var1 var2 cor statistic p method
## <chr> <chr> <dbl> <dbl> <dbl> <chr>
## 1 percdead Dim.2 -0.011 4104. 0.955 SpearmanGraphs
Tests
percdead_long<- percdead_glm_data %>% dplyr::select(percdead,spawn_nh4n,spawn_no3n,forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain) %>% pivot_longer(.,cols=c(spawn_nh4n,spawn_no3n,forest_area,arable_area,grassland_area,freshwater_area,urban_area,other_area,tmax,tmin,af,rain),names_to="variable",values_to="traitval")
percdead_corr<-percdead_long %>% group_by(variable) %>% cor_test(traitval,percdead,method="spearman")
percdead_corr<-mutate(percdead_corr,p.adj=p.adjust(p,"BH"))
percdead_corr$Variable<-label_df$new[match(percdead_corr$variable,label_df$old)]
percdead_corr%>% dplyr::select(Variable,cor,statistic,p,p.adj) %>% dplyr::arrange(p.adj)## # A tibble: 12 × 5
## Variable cor statistic p p.adj
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Arable Area -0.53 6193. 0.00343 0.0209
## 2 Maximum Temperature -0.52 6189. 0.00349 0.0209
## 3 Grassland Area 0.49 2077. 0.00719 0.0288
## 4 Urban Area -0.37 5571. 0.0469 0.141
## 5 Spawn Nitrate -0.35 5481. 0.0627 0.150
## 6 Minimum Temperature -0.31 5323. 0.101 0.202
## 7 Forest Area -0.28 5192. 0.143 0.245
## 8 Air Frost 0.26 3020. 0.18 0.27
## 9 Other Area 0.18 3346. 0.361 0.481
## 10 Freshwater Area 0.16 3419. 0.414 0.497
## 11 Spawn Ammonium -0.069 4342. 0.72 0.785
## 12 Rain 0.0099 4020. 0.959 0.959Graphs
#Arable area
percentdead_plot2 <- ggplot(percdead_df, aes(arable_area/1000000, percdead)) +
xlab(bquote("Arable area " (km^2))) +
ylab("Percent Dead Spawn") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location)) +
theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + theme(legend.position = "top") + annotate(geom="text",x=7,y=75,label="cor = -0.53 \np.adj=0.02",size=6)
percentdead_plot2## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
#TMax
percentdead_plot3 <- ggplot(percdead_df, aes(tmax, percdead)) +
xlab("Mean Maximum Temperature") +
ylab("Percent Spawn \n Dead (Logits)") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location)) +
theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + theme(legend.position = "top") + annotate(geom="text",x=7,y=75,label="cor = -0.52 \np.adj=0.02",size=6)
percentdead_plot3 #+ facet_wrap(.~location)## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
#Grassland Area
percentdead_plot4 <- ggplot(percdead_df, aes(grassland_area/1000000, percdead)) +
xlab(bquote("Grassland area " (km^2))) +
ylab("Percent Dead Spawn") +
geom_smooth(method = "lm") + geom_point(shape=21,size=5,aes(fill=location)) +
theme_bw(base_size = 15) + labs(fill="Location") + scale_fill_brewer(palette = "Set2") + theme(legend.position = "top")+ annotate(geom="text",x=10,y=75,label="cor = 0.49 \np.adj=0.03",size=6)
percentdead_plot4## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
####### SAVE GRID
# percdead_grid1<-plot_grid(percdead_rotationplot,percentdead_plot2,ncol=2,labels="AUTO",label_size=25)
# percdead_grid2<-plot_grid(percdead_grid1,percdead_site_plot2,nrow=2,labels=c("","C"),label_size = 25)
#
# ggsave2('Fig_5_Percent Dead Combined Plot.pdf',percdead_grid2,width=25,height=25,units="cm")Individual Predictors
#Save PCA Plot
ggsave2('Percent Dead PCA PLot.pdf',percdead_rotationplot,width=10,height=10,units="cm")
#Individual Predictors
percdead_combinedplot<-plot_grid(percentdead_plot2,percentdead_plot3,percentdead_plot4,labels=c("C","D","E"),nrow=2,label_size = 25)## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).#Stich On Site Specific Temp Trends
percdead_combinedplot2<-plot_grid(percdeaad_siteplot1_ann,percdead_site_plot2,ncol=2,labels = c("A","B"),label_size = 25)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'#Stich On Bayesian Model
percdead_combinedplot3<-plot_grid(percdead_combinedplot2,percdead_combinedplot,nrow=2,rel_heights = c(1.5,2))
#Save
ggsave2('Fig_4_PercentDead Grid Plot.pdf',percdead_combinedplot3,width=40,height=45,units="cm")